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1.  Introduction.  The  development  of  computers  thirty  years  ago  made 
it  practical  to  calculate  finite  difference  approximations  of  elliptic 
partial  differential  equations.  For  these  calculations  the  solution  of 
a  linear  system  AU  =  F,  which  is  the  finite  difference  representation 
of  the  differential  equation,  is  fundamental.  Hardware  characteristics 
of  early  computers,  particularly  memory  limitations,  spurred  the 
development  of  direct  iterative  methods  for  these  linear  systems.  In 
direct  iterative  schemes  the  matrix  A  splits  into  a  difference  A  =  M-N, 
and  one  generates  a  sequence  {U^}  according  to  MU^  =  NU^V  ^+F. 
Convergence  of  the  method  is  governed  by  the  spectral  radius  p  of  M  *N: 
{U^}  converges  to  the  solution  if  p  <  1,  and  smaller  p  implies  faster 
convergence. 

The  first  iterative  methods  were  point  methods  —  in  any  step  of 
the  iteration  they  solved  for  one  component  of  the  unknown  solution 
vector  at  a  time.  Intuition  suggests  that  iterative  algorithms  that 
solve  for  several  points  at  once  will  converge  more  rapidly  than  point 
algorithms.  The  Gaussian  elimination  algorithm  is  seen  in  this  light 
to  converge  in  one  step.  Frankel  l 14],  Young  [34],  Arms,  Gates,  and 
Zondek  [1],  and  Varga  [32],  using  the  algebraic  structure  of  the  linear 
systems,  and  Parter  [22],  [23],  by  exploiting  the  nature  of  the  systems 
as  finite  difference  approximations  to  elliptic  partial  differential 
equations,  determined  the  convergence  rates  of  point  and  block 
iterative  methods.  The  results  confirmed  that  iterative  methods  on 
blocks  comprising  several  lines  of  unknowns  indeed  converged  faster 
than  point  methods.  Much  of  the  work  up  to  1961  is  collected  in  [33], 

The  usual  finite  difference  approximations  are  accurate  to  second 
order  in  the  spatial  mesh  size  h.  In  the  middle  1960s  attention  turned 
to  higher  order  approximation  methods  —  finite  element  and  other 
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projection  methods,  which  are  still  the  subject  of  intensive  study 
(136),  [7],  [2],  [30],  [5],  [10]).  Because  of  their  treatment  of 
boundary  conditions,  these  methods  are  formally  easier  to  obtain  than 
higher  order  finite  difference  approximations,  and  for  a  given  accuracy 
their  corresponding  linear  system  of  equations  is  smaller  than  the 
finite  difference  system.  Hence  interest  in  direct  factorization 
methods  for  linear  systems  grew,  and  continues  today;  see  [27],  [28], 
[15],  and  [16]. 

At  about  the  same  time  it  was  seen  that  their  regular  structure 
made  separable  finite  difference  elliptic  systems  amenable  to  special 
fast  direct  factorization  methods  ([18],  [9],  [12],  [31]).  For  a 
limited  class  of  nice  elliptic  problems,  then,  it  became  practical  to 
compensate  for  the  second  order  accuracy  of  the  usual  finite  difference 
approximation  by  taking  a  sufficiently  small  h  and  exploiting  the 
regular  structure  of  the  linear  system. 

But  not  every  problem  is  nice.  Moreover,  within  the  past  few 
years  a  growing  desire  to  solve  three-dimensional  problems,  together 
with  the  development  of  novel  computer  architectures  —  array 
processors,  vector  machines,  and  multiprocessors  —  has  rekindled 
interest  in  block  iterative  methods  for  elliptic  systems.  The  effects 
of  special  architectures  are  considered  in  [29],  [17],  and  [19],  while 
an  analysis  of  the  convergence  rates  of  iterative  methods  for  fairly 
general  elliptic  problems  already  appears  in  [23]. 

But  not  every  analysis  is  nice,  and  that  of  [23],  partly  because 
of  its  generality,  is  somewhat  opaque.  A  relatively  direct  discussion 
of  the  basic  ideas  is  given  in  [3]  for  the  Poisson  problem  in  a  square. 
That  presentation  uses  the  strong  estimates  of  Nitsche  and  Nitsche  [21] 
and  of  Brandt  [8]. 


Our  purpose  here  is  to  reexamine  the  convergence  rates  of 

A 

iterative  block  methods  for  elliptic  difference  equations.  A  feature 
of  the  present  analysis  is  that  we  avoid  the  estimates  of  {21]  and  (8). 
For  the  Poisson  problem  in  two  or  three  dimensions  this  is  of  little 
moment.  But  the  Nitsche  estimates  have  never  been  extended  to  general 
regions,  and  must  fail  in  dimensions  greater  than  three.  In  contrast, 
we  will  show  that  our  new  approach  is  easily  extended  to  general 
domains,  to  any  number  of  dimensions,  and  to  general  elliptic 
difference  equations. 

In  addition,  we  can  deal  with  certain  kinds  of  singularly 

perturbed  elliptic  difference  equations.  Such  equations  can  arise  when 

solving  parabolic  problems  by  discrete  time  methods.  For  instance,  let 
2  2  2 

A  :=  3  /3xi  be  the  two-dimensional  Laplacian;  the  backward  Euler 

method  for  the  parabolic  operator  (Cp3/3t)-A  leads,  at  each  time  slice 
tQ,  to  an  elliptic  operator 

(1.1)  c/x  -  A,  T  :=  tQ  -  tQ_j. 

Let  A^  be  a  finite  difference  approximation  to  A  on  a  spatial  mesh  of 

size  h;  we  get  a  matrix  A  representing  the  elliptic  difference 
2  2  2  -0 

operator  ch  /x-h  A^.  If  ch  /x  =  ch  ,  then  A  corresponds  to 

(1.2)  ch“  -  h2A^. 

We  distinguish  four  cases.  Analysis  of  the  first,  in  which  a  <  0, 
is  easy:  p  =  0(h"°),  and  iterative  methods  converge  very  rapidly.  In 
the  second,  a  =  0,  and  (1.2)  is  a  singularly  perturbed  operator.  We 
have  studied  this  operator  in  [26],  where  it  arose  from  plane  iterative 
methods  for  the  Poisson  problem  in  the  unit  cube;  the  attack  there, 
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though  related  to  some  of  the  ideas  of  this  report,  seems  to  be 
particular  to  the  model  operator  (1.1)  and  rectangular  domains. 

In  this  paper  we  restrict  our  attention  to  the  third  and  fourth 
cases,  wherein  0  <  a  £  2.  If  a  -  2,  then  (1.2)  is  a  regular  elliptic 
difference  operator,  to  which  both  the  earlier  and  our  new  analyses 
apply.  When  0  <  a  <  2,  (1.2)  is  again  a  singularly  perturbed  operator; 
but  it  too  can  be  handled  with  our  present  methods,  unlike  the  instance 
or  =  0.  To  justify  considering  this  case,  we  point  out  that  a  =  1  for 
the  optimal  choice  of  x  in  the  Crank-Nicolson  method  for  parabolic 
problems. 

We  begin  in  section  2  with  a  description  of  the  model  elliptic  and 
parabolic  problems  in  the  two-dimensional  unit  square.  It  is  worth 
remarking  that  our  model  problems  need  not  be  self-adjoint.  Section  3 
is  devoted  to  proving  the  convergence  rates  of  iterative  schemes 
satisfying  certain  basic  assumptions. 

In  section  4  we  describe  block  structures  of  particular  interest 
—  k-line  and  k*k  blocks  —  and  the  usual  iterative  schemes:  Jacobi, 
Gauss-Seidel ,  and  successive  overrelaxation.  In  these  schemes  A  splits 
into  a  difference  A  =  M-N.  The  key  to  our  analysis  is  that  it  suffices 
to  consider  only  the  block  Jacobi  scheme,  for  which  N  is  essentially  a 

a* 

sum  of  one-dimensional  weak  multiplication  operators  N.  We  demonstrate 

A* 

this  decomposition  of  N  in  section  4,  and  discuss  the  action  of  N  in 
section  5.  In  section  6  we  use  the  theory  of  section  3  and  properties 
of  N  to  derive  the  convergence  rates  of  the  block  iterative  methods  of 
section  4. 

Next  we  take  up  more  general  problems:  other  operators  in  section 
7,  and  other  domains  in  section  8.  We  conclude  in  section  9  with  some 
comments  about  the  general  applicability  of  our  method  of  analysis. 
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2.  The  model  problems.  The  basic  ideas  are  clearest  in  this  simple 
setting.  We  construct  finite  difference  approximations  of  the  partial 
differential  operators 


Lu  :=  -Au  +  du  +  eu  , 
x  y 

(2.1) 

£u  :=  cut  +  Lu 
on  the  open  unit  square 

0  :=  {(x,y)  €  R2  :  0  <  x,  y  <  1} 
in  the  usual  way.  Impose  on  0  a  mesh  with  uniform  spacing 
(2.2)  h  :=  1/(P  +  1) 


and  let  (x^y^)  :=  (ih,jh).  Define  the  set  of  interior  mesh  points  0^ 
and  the  discrete  boundary  30^  by 


(2.3) 


:=  {(xi,yj)  :  1  £  i,  j  £  P}, 

30.  :=  {(x.,y.)  :  i  =  0  or  =  P+1,  or  j  =  0  or  =  P+1 } 
n  i  j 


A  mesh  vector  U={U.  .  :  0  £  i,  j  $  P+1}  is  a  function  defined  on  the 

1 1 J 

entire  discrete  mesh  0^  :=  0^  u  30^. 

The  discrete  Laplace  operator  is  defined  at  points  in  0^  by 

'V'i.j  !=  '  2”i,j *  W"2 

(2.4) 

+  '  2Ui.j  * 

We  suppose  that  c,  d,  and  e  are  smooth  functions  on  0  and  that 


(2.5)  c(x,y)  £  cn  >  0  on  0. 
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The  discrete  operators  that  arise  in  approximating  (2.1)  are  then 


(2.6) 


[L.  Uj .  .  :=  1-A.Uj.  .  +  d.  .(U.A1  .  -  U.  .  .)/(2h) 
h  Ji,j  Ti  Ji,j  i,j  l+l, j  i-l,j"v 


+  e.  .(U.  -  U.  .  . )/(2h) 

i.J  i,j+l  i.j-l"'  . 


(2.7)  U,U],  .  :=  (c.  ,/x)U.  .  +  [L.U] .  ., 

h  i,J  i.J  i,J  b  i>j’ 

where  X  >  0  is  given  and,  for  instance,  c.  .  :=  c(x.,y.). 

i » J  i  3 

Note  that,  although  the  mesh  vector  U  is  defined  on  the 
vectors  A^U,  L^U,  and  £^U  are  defined  only  at  the  interior  mesh  points. 
As  usual,  the  forward  difference  operators  are  given  by 


Vi.J  <Ui+l,j  *  Ui,j>/h  (0  *  i  s  p,  1  S  JS  P), 

(2.8) 

VyUi,j  **=  (Ui,j+l  “  ui,j)/h  (^i^P.  OSJSP). 

Given  mesh  vectors  F  and  G,  the  model  elliptic  problem  is  to  find 
a  mesh  vector  U  satisfying 

(2.9)  LjU  =  F  in  U  =  G  on  3^ 

and  the  model  parabolic  problem  requires  U  to  solve 

(2.10)  ijU  =  F  in  fib,  U  =  G  on  3^. 


After  choosing  an  ordering  of  the  mesh  points  (x.,y.)  -«  or, 

i  J 

equivalently,  of  the  components  of  U  —  we  let  A  be  the  matrix 
2  2 

representing  h  or  h  2^.  As  indicated  in  (2.4),  A^,  L^,  and  map 

2  2 
vectors  with  P  +4P  components  into  vectors  with  P  components.  Hence  A 

2 

is  a  matrix  of  order  P  ;  the  known  boundary  values  G  are  put  on  the 
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right  hand  sides  of  the  difference  equations  (2.9)  and  (2.10).  In 
either  case  we  arrive  at  a  linear  system 


(2.11)  AU  =  F 

2  ~ 

of  order  P  ,  where  F  indicates  the  result  of  ordering  the  components  of 
2 

h  F  and  of  including  the  G  terms. 

2 

Every  vector  U  with  P  components  may  be  viewed  as  a  mesh  vector 
on  that  also  satisfies 

(2.12)  U  =  0  on  3^. 

Henceforth  we  assume  every  mesh  vector  U  satisfies  (2.12) . 

An  iterative  method  for  solving  (2.11)  is  determined  by  a 
splitting 

(2.13)  A  =  M  -  N. 

Rewrite  (2.11)  as 

MU  =  NU  +  F. 

After  choosing  a  first  guess  U^,  we  obtain  a  sequence  {U^}  from 

(2.14)  MU(V)  =  HU(V_1)  +  F. 

It  is  well  known  that  when  A  is  nonsingular  the  iterates  {U^} 
converge  to  the  unique  solution  of  (2.11)  independently  of  if  and 

only  if  the  spectral  radius 

p  :=  max  {|A|  :  det(XM-N)  =  0) 

of  M~^N  satisfies  p  <  1.  So  the  first  thing  we  require  of  a  splitting 


is  that  p  <  1.  Evidently  the  iterates  {IT  of  (2.14)  converge  more 
rapidly  for  smaller  p.  Hence  our  task  is  ‘to  determine  the  asymptotic 
behavior  of  p  as  h  -*■  0. 

For  future  reference  we  note  that  corresponding  to  every  X  for 
which  det(XM-N)  =  0  there  is  a  vector  V  t  0  satisfying  XMV  =  NV.  We 
also  record  two  lemmas  regarding  V^,  V  ,  and  A^.  Let  X  and  Y  be  mesh 
vectors;  define  an  inner  product  and  associated  norm 


(X,Y)  :=  2.  .  X.  .Y.  IX  B.  :=  (X.X)1^. 
x,j  i,j  i,j’  h 


An  operator  B  on  mesh  vectors  is  normed  in  the  customary  way  by 


IB I.  :=  sup  {lBX«h  :  »Xlh  =  l}. 

As  usual,  | d | ^  denotes  the  sup  norm  of  d  over  5. 


Lemma  2.1.  If  U  is  a  mesh  vector  satisfying  (2.12),  then 


(VU,VU)  +  (VU.VU)  =  (~A^U,U) 


Proof.  Summation  by  parts;  see  [11]  or  [20],  D 


Lemma  2.2.  If  U  is  a  mesh  vector  satisfying  (2.12),  then 


(|VxU|,|u|)  +  (|V  U|,|U|)  S  IUIh[2(-\U,U)]‘ 


Proof.  By  the  Schwarz  inequality. 


(ivui.iun  +  (|vyu|,|u|)  $  iuih[ivuih  +  ivyuih], 

2  2 

But  the  inequality  2ab  £  a  +b  and  Lemma  2.1  show  that 


[,yj.h  +  lVyU8h]2  S  2|8V"h2  +  ,VyU"h2]  =  2(-AjjU,U) .  □ 
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3.  A  general  approach.  To  begin  the  analysis  of  the  splitting  (2.13) 


we  make  four  assumptions. 

Al.  p  <  1,  so  the  iterative  method  (2.14)  is  convergent. 

A2.  p  is  an  eigenvalue  of  M  there  is  a  mesh  vector  U  *  0  such 

that  pMU  =  NU. 

A3.  There  is  a  positive  constant  Nq,  independent  of  h,  such  that 
INIhSN0. 

A4.  There  are  a  smooth  function  q  and  constant  with 
q(x,y)  ^  q0  >  o  on  q 

and  a  constant  D  >  0,  independent  of  h,  so  that  whenever  U  and  V 
are  mesh  vectors  satisfying  (2.12)  we  have 

(NU,V)  =  (qU,V)  +  E, 

where 


|E|  $  hDKlV  U|  +  |V  U|,|V|)  +  (|U|,|V  V|  +  |V  V|)  +  (|U|,|V|)] 

*  y  *  y 

+  h2D[(-AjU,U)  +  (-AjV/V)]. 

Assumptions  Al  -  A3  are  in  effect  more  or  less  common;  this  will 
become  clear  in  section  6.  Our  main  new  concept  is  A4.  As  might  be 
expected,  verification  of  A4  and  the  determination  of  q  are  the 
important  technical  steps  when  applying  our  analysis  to  any  particular 
splitting.  But  we  shall  see  that  these  steps  are  not  difficult. 

10 
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When  a  splitting  (2.13)  satisfies  these  assumptions,  the 
asymptotic  behavior  of  p  as  a  function  of  h  is  readily  discovered.  We 
begin  with  the  elliptic  case. 

2 

Theorem  3.1.  Let  A  correspond  to  h  L^.  Suppose  the  splitting 
(2.13)  satisfies  A1  -  A4.  Let  be  the  smallest  eigenvalue  of  the 
problem 

(3.1)  Lv  =  Aqv  in  Q,  v  =  0  on  3fl. 

Then 

(3.2)  p  =  1  -  AQh2  +  o(h2). 

Proof.  Let  U  be  the  eigenvector  associated  with  p  in  A2,  so  that 


pMU  =  NO. 

Subtract  pNU  from  both  sides  and  use  (2.13)  to  see  that 

(3.3)  AU  =  ((1  -  p)/p)NU. 

By  Al, 

(3.4)  jj  :=  (1  -  p)/(ph2) 

2 

is  positive.  Because  A  represents  h  L^,  (3.3)  corresponds  to 

(3.5)  LhU  =  pNU  in  U  =  0  on  3nfa. 

Indeed,  whenever  A  t  0  satisfies 

(3.6)  AMX  =  NX 

for  some  nonzero  X,  then 
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p  =  p(A)  :=  (1  -  M/(^2) 


is  an  eigenvalue  of  (3.5).  Conversely,  if  p  is  an  eigenvalue  of  (3.5) 

2 

and  1+ph  t  0,  then 

\  =  X(M)  :=  1/(1  +  yh2) 

is  an  eigenvalue  of  (3.6). 

For  fixed  h,  let  p  be  an  eigenvalue  of  (3.5)  minimal  in  magnitude. 

The  basic  result  of  [24]  shows  that  p  ■+  Aq  as  h  ■+  0  —  that  is,  p  = 

.  2 

Aq+o(1).  It  follows  by  positivity  of  Aq  that  Re(l+ph  )  >  0  for  small 
h,  whence  A  :=  l/(l+ph  )  is  a  well  defined  eigenvalue  of  (3.6).  Hence 

p  Z  |A|  =  1/ | 1  +  ph2 |  =  1  -  [Aq  +  o(l)]h2. 

But  p  given  by  (3.4)  is  an  eigenvalue  of  (3.5)  by  construction,  and  so 
(l-p)/(ph2)  Z  Ipl  =  Aq+o ( 1 ) ,  by  the  minimality  of  p.  We  deduce  that 

p  i  1/(1  +  [AQ  +  o(l)]h2)  =  1  -  [Aq  +  o(l))h2. 

Comparison  of  this  and  the  previous  inequality  proves  (3.2).  □ 

Parabolic  equations  lead  to  discrete  singular  perturbation 

eigenvalue  problems,  so  in  the  general  nonself-adjoint  case  we  can 

establish  only  an  inequality  analogous  to  (3.2).  We  arrive  as  before 

2 

at  (3.3),  where  A  represents  h  2^;  hence 
(3.7)  h2(c/x  +  Lfa)U  =  ((1  -  p)/p)NU. 

We  make  a  basic  assumption  about  the  ratio  of  the  time  step  T  to  the 
spatial  mesh  size. 
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PI.  There  are  constants  >  0  and  0  <  a  <  2  such  that  h2/x  =  c^h01. 

k 

Now  define 

(3.8)  M  :=  (1  -  p)/(pha); 

we  deduce  from  (3.7),  (3.8),  and  PI  that  in  the  parabolic  case  (3.3) 
corresponds  to 

(3.9)  cU  +  h2-0^  =  pNU  in  f^,  U  =  0  on  3^, 
where 

(3.10)  c(x,y)  :=  c1c(x,y). 

2 

Theorem  3.2.  Let  A  correspond  to  h  2^.  Suppose  PI  holds  and  the 
splitting  (2.13)  satisfies  A1  -  A4.  Let 

Aj  :=  min  {c(x,y)/q(x,y)  :  (x,y)  €  fi}. 

Then 

(3.11)  pSl-  Ajh0  +  o(h0). 

Proof.  Because  p  is  positive,  (3.11)  is  equivalent  to  (l-p)/(pha) 
=  p  £  Aj+o(l).  Suppose  this  inequality  is  false.  We  may  then  assume 

(3.12)  0  $  p  S  2Aj . 

Let  U  be  the  eigenvector  of  (3.9)  associated  with  p.  Normalize  1131^  to 
be  1.  By  A4, 

(3.13)  p(NU,U)  =  p(qU,U)  +  E1? 


where,  using  Lemma  2.2, 


|Ej|  %  2|JhD[2(-AhU,U)]1/2  +  2ph2D(-AhU,U)  +  phD. 

Use  (3.12)  and  the  inequality  2ab  %  a20  2+b202  to  get 

|Ej|  <  16A12D2h“  +  B/2  +  AAjDh^  +  2A^Dh, 
where  we  have  defined 

B  :=  h2"a(-AhU,U). 

Lemma  2.1  shows  that  B  >  0.  It  follows  from  (3.9)  and  (3.13)  that 
(cU,U)  +  B  =  M(qU,U)  +  Ej  +  e2, 

with 

|E2I  5  h2"“K(2(-AhU,U)]1/2  S  2h2'aK20"2  +  B02 

and  K  :=  | d | oo+ 1 e 1 ^ .  Choose  0  so  small  that  the  coefficients  of  B  in 
the  estimates  of  Ej  and  E2  sum  to  less  than  1  for  small  h.  Then 

(cU,U)  3  p(qU,U)  +  2h2'aK20’2  +  +  2AjDh. 

The  theorem  follows  at  once.  □ 

When  the  splitting  is  self-adjoint  —  a  frequent  occurrence  —  we 
can  use  the  variational  principle  to  establish  equality  in  (3.11). 

Theorem  3.3.  Under  the  assumptions  of  Theorem  3.2,  suppose  also 
that  we  have 

SI.  A  and  M  are  Hermitian  and  positive  definite. 

Then 

(3.14)  p  =  1  -  Ajh0  +  o(h°). 
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r 


4  - 

Proof.  Fix  €  >  0  and  choose  v  (x,y)  e  C  (Q)  to  vanish  on  3Q  and 

V 


to  satisfy 


(3.15) 


Jn  cv£2  dx  dy 
;fl  qve2  dx  dy 


£  Aj  +  C- 


Now  A2  and  SI  imply  that  p  =  sup  {(NX,X)/(MX,X)  :  X  t  0}.  Choosing  X 
as  the  mesh  vector  V  determined  by  point  evaluation  of  v  yields 

(3.16)  p  £  (NV£,V£)/(MV£,V£)  =  (NV£,V£)/((AV£,V£)  +  (NV£,V£)J. 

Observe  that 


«W  *  *  h2'“(1hVV>' 

It  follows  from  the  smoothness  of  v  that 

e 

h2(AW  =  cv£2  dx  dy  +  0(h2’a)l; 

moreover,  by  A4 

h2 (NV£,V£)  =  qv£2  dx  dy  +  o(l). 

Combining  these  equalities  with  (3.15)  and  (3.16)  yields 
P  il  -  (Aj  +  e)h“  +  o(ha), 

which  together  with  (3.11)  establishes  the  theorem.  Q 


Note  that  hypothesis  SI  requires  d  =  e  =  0  for  the  operators  (2.1) 
of  the  model  problems. 
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4.  Some  block  iterative  methods.  We  take  up  now  a  description  of 
specific  block  iterative  methods  corresponding  to  (2.13).  The  block 
structure  of  an  iterative  scheme  for  the  linear  system 

AX  =  Y, 


where  A  is  an  n><n  matrix,  is  completely  determined  by  a  block  partition 
of  the  n-vectors.  Suppose  every  n-vector  X  is  decomposed  into 
subvectors 

x  =  (xrx2,  ..  ,xr)fc 

and  each  X^  is  itself  an  n^ -vector.  This  partition  of  X  induces  a 

block  partition  A  =  [A.  .]  in  which  each  A.  .  is  an  n.xn.  matrix.  The 

1 » J  1  *J  1  J 

corresponding  block  Jacobi  iterative  scheme  is 

(4.1)  A.  .X.(V)  =  -  I  A.  X  (V-1)  +  Y.. 

1,1  1  S*1  1,6  S  1 

In  terms  of  (2.14),  M  is  the  block  diagonal  matrix  M  =  diag[A.  .].  The 

1,1 

corresponding  Gauss-Seidel  scheme  is 

(4.2)  A.  ,X.(v)  =  -  Z  .  A.  X  -  Z  .  A.  X  (v_1)  +  Y., 

1,1  1  S<1  1,S  S  S>1  l,S  S  l’ 

while  the  successive  overrelaxation  (SOR)  method  with  relaxation 
parameter  in  is 


(4.3) 


A.  .X.(W)  =  -  uuZ  ..  A.  X 

1,1  1  S<1  1,S  s 


wZ  A.  X  -  uiZ  v.  A.  X 


(v-1) 


+  U)Y.  +  (1  -  w)A.  .X. 


s>i  i,s  s 

(v-1) 


We  are  interested  in  specific  block  structures  that  arise  in  a 


natural  geometric  way.  Recall  that  a  mesh  vector  U  is  defined  on  the 
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rectangular  set  of  mesh  points  f}^.  We  will  decompose  U  into  blocks  of 
components  corresponding  to  lines  or  subsquares  of  mesh  points. 
Formally,  let  k  be  a  fixed  integer  factor  of  P,  so  that 

(4.4)  P  =  kQ  for  some  integer  Q. 


In  the  k-line  block  structure  (see  [22]  or  [23]  for  a  detailed 
description),  each  block  of  U  comprises  the  unknowns  IL  ^  associated 
with  the  points  on  k  consecutive  horizontal  (or  vertical)  grid  lines. 
Index  the  blocks  by  s;  we  have 

(4.5)  0S  :=  ••  1  S  i  4  P.  1  S  j  S  «. 

The  k*k  block  structure  is  described  in  [3],  [25],  and  [26].  Each 
block  comprises  the  unknowns  associated  with  a  kxk  square  of  mesh 
points.  We  distinguish  these  blocks  with  a  double  index  (r,s): 


(4.6) 


Ur,s  *“  ^Uk(r-l)+i,k(s-l)+j 


1  S  i,  j  £  k}. 


To  write  down  the  matrices  A,  M,  and  N  of  the  Jacobi  iterative 
method  for  each  of  these  block  structures  is  straightforward  but 
tiresome.  We  shall  give  a  unified  analysis  of  the  Jacobi  method  for 
these  structures.  But  for  illustrative  purposes  we  first  sketch  a 
development  of  the  (horizontal)  k-line  scheme  for  the  elliptic  problem 
(2.9). 

For  1  $  0  S  P  define  the  P*P  matrices 


(4.7)  Sc  :=  diag|l*h«ij0/2]  (1  S  i  S  P) 
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The  notation  indicates  that  is  tridiagonal  while  SQ  and  Ta  are 
diagonal.  For  example, 


'Vi  .r< 


0  if  |i  -  jl  >  1 

•l-hdi>0/2  if  j  =  i  -  1 


4  if  j  =  i 

|^-l+hdi>0/2  if  j  =  i  +  1. 


With  this  ordering  of  the  mesh  points  into  horizontal  lines,  A  is  the 
2  2 

P  xp  block  tridiagonal  matrix 


(4.8)  A  =  f-S0,  Dff,  -T0]  (UffSP), 


Now  collect  the  lines  of  unknowns  k  at  a  time.  For  1  ^  s  S  Q  let 
be  the  kPxkP  block  tridiagonal  matrix 

Ms  :=  [-Sk(s-l)+CT’  Dk(s-l)+a*  “Tk(s-l)+o^  (1  =  a  =  k)’ 


and  define  the  kPxkP  block  matrices 


R  :  = 


Tks  0 


W  :  = 


0  Sk(s-1)+1 

0  0 


Observe  that  A  is  then  the  block  tridiagonal  matrix 


A  =  [-Ws,  Ms,  -RsJ  (1  <  s  i  Q). 


In  the  k-line  Jacobi  scheme,  A  splits  into  the  block  matrices 


M:=diag(Ms],  R  :*  lWgl  0,  »g), 


We  now  seek  a  simple  quantitative  description  of  N  for  both  the 


k-line  and  the  kxk  block  partitions  when  k  £  2.  If  B  and  C  are 


matrices,  we  mean  by  B  =  C+0(h)  that  there  is  some  constant  K  so  that 

a 

| (BX,Y)  -  (CX,Y)|  £  Kh| (X,Y) |  for  every  X  and  Y. 

Because  and  T0  are  0(h)  perturbations  of  the  PxP  identity  matrix, 
let  us  for  the  moment  ignore  the  small  terms.  We  define  a 
one-dimensional  operator  N  on  vectors  <J>  :=  (<|>^,  (j)^,  ..  ,$p)t  as 
follows: 


(4.9) 


[N4>] 


ks+o 


:= 


*ks+l 


'ks 


1  S  s  i  Q-l,  a  =  0 

1  S  s  %  Q-l,  o  =  1, 


0  for  any  other  subscript  j. 


N  is  a  weak  multiplication  operator,  as  we  shall  see  in  the  next 
section.  Now  let  Nx  be  that  operator  on  mesh  vectors  U  that  acts  on  U 
only  in  the  x-direction,  and  in  that  direction  acts  as  N.  Define  in 
a  similar  way.  For  instance,  with  1  £  i  ^  P  we  have 


lyii.i.*, :=<; 


(4.10) 


^i,ks+l 


Ui,ks 


lNyu]i,i  := 


1  i  s  i  Q-l,  a  =  0 

1  S  s  S  Q-l,  o  =  1, 

for  any  other  subscript  j . 


Observe  for  each  block  structure  that  the  Jacobi  splitting  (4.1) 
yields  the  same  N  for  both  the  elliptic  and  parabolic  operators  (2.6) 
and  (2.7).  This  is  so  because  the  matrix  representing  the  operator 
4^-1^  is  a  diagonal  matrix.  A  straightforward  computation  proves  the 
next  theorem,  which  summarizes  the  essential  nature  of  N. 
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Theorem  4.1.  Let  k  i  2.  In  the  k-line  Jacobi  scheme  (4.1)/(4.5), 


(4.11)  N  =  Ny  +  0(h) , 

and  for  the  k*k  block  Jacobi  scheme  (4.1)/(4.6)  N  is  given  by 

(4.12)  N  =  N  +  N  +  0(h).  D 

x  y 


5.  The  operator  N.  We  now  show  that  NU  converges  weakly  to  (2/k)U,  so 
that  N  is  a  weak  multiplication  operator.  In  this  section  U  and  V  are 
real  vectors  with  P  components.  For  each  such  vector  X  it  is  useful  to 
define  Xq  :=  0.  It  is  clear  from  (4.9)  that  N  samples  U  twice  in  each 
block  of  k  points  {U.  ^  :  0  £  a  %  k-l},  where  0  £  s  %  Q-l  --  except  in 
the  first  and  last  blocks.  Roughly,  but  perhaps  vividly,  N  sees  U 
about  2/k  of  the  time;  precisely,  from  (4.9)  we  have 

(5.D  cira.v)  =  1%1  (uksvkstl  ♦  <Jkstlvk5). 

If  U  and  V  arise  from  the  evaluation  of  smooth  functions  u(x)  and 
v(x)  on  the  points  {x^  :=  ih  :  0  ^  i  i  P+1}  ,  then 


U.  ..  £  U,  and  V,.  =  V, 


ks+j 


ks 


ks+j  ~  ks+1 


(0 


k-l), 


whence 


Vte.l=Wb*i  (0  £  j  5  k-l). 

Summing  this  approximate  equality  over  j  and  dividing  by  k  gives 

VWl  5  'WW 


which  looks  like  a  Riemann  sum  over  the  interval  [x^  , 
(1/hk)/  u(x)v(x)  dx.  Consequently, 


xks+k^ 


for 


(NU,V)  S  (2/hk)  /[0)1]  u(x)v(x)  dx  S  (2/k)(U,V). 

Now  we  make  this  argument  precise.  Let  V  be  the  forward 
difference  operator,  as  in  (2.8).  Fix  j  for  the  moment.  Obviously 

Uks  =  Uks+j  "  ^0=0  ^ks+o'  Vks+1  =  Vks+j  '  ^0=1  Wks+o 
(as  usual,  a  vacuous  sum  has  value  0).  Hence 


(5.2) 


°ksVks+l  *  Wk.+J  +  »2‘Zi£j  TOkS4a)(Zi£l  OTks4„> 
•  W'Wo  -  “£2  Vks+jVUks+a' 


Replace  Ujts+j  and  the  last  two  terms  of  (5.2),  using  the 

identity 


X,  ,  *-  X,  .  .  -  hljli  VX.  .  . 

ks+o  ks+j  n=a  ks+n 


This  substitution  gives 


VW  =  °k.4jvkS4j  +  h  Go,j<»Gi,j<v> 

(5.3)  -  h£|  Ukst0Wks+o  -  kljlj  Vks+0VUksM 

-  “2C!  G<,,.j™Wks«r  -  ih£  Ga,j<V>'rai 'ks«. 
where  for  0  &  a  %  j-1  £  k-1  we  define 

Ga,j®  -  z£  rak,4„.  «*■>>  :=  £2  "W  «  lGc,j<«' 
Sum  (5.3)  over  0  5  j  £  k-1  and  divide  by  k  to  get 

<5-‘>  Vks.l  *  «'k>ZS2  Gks+jEks+j  *  <1/k)Es' 
with 
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(5.5) 


|E  |  i  hk[Zk'];  |U,  ..IIW,  A.|  +  Zk_*  |V,  ^ .  1 1 VU,  ^.|] 
s'  j=0  ks+j  ks+j '  j=0  ks+j  1  ks+j 1 1 

* 

+  3h2kG(U,s)G(V,s). 

2  2 

By  the  Schwarz  inequality  and  the  inequality  2ab  %  a  +b  , 

G(U,s)G(V,s)  S  (k/2){(ZklJ  IWks+j|2)  +  (ZkZj  |Wks+j|2)]. 

Estimate  the  last  term  of  (5.5)  in  this  way,  and  sum  (5. A)  over  s  to 
deduce  that 


l%l  UksVks*l  =  0/k)(U,V)  ♦  E/2, 

where 


(5.7) 


|E|  *  2h[(|U|,|W|)  +  ( |VU| ,  |V|  )] 
+  3h2k[  (VUjVU)  +  (W,W)}. 


Comparison  of  (5.1)  to  (5.6)  shows  that  we  can  exploit  the 
symmetry  of  this  argument  in  U  and  V  to  prove  the  following  theorem, 
which  quantitatively  describes  N. 


Theorem  5.1.  Let  N  be  given  by  (4.9).  For  P-vectors  U  and  V, 
(5.8)  (MJ,V)  =  (2/k) (U,V)  +  E, 

w  f 

and  E  is  estimated  by  (5.7).  □ 


6.  Rates  of  convergence.  In  this  section  we  take  up  the  problem  of 
determining  the  convergence  rates  of  the  iterative  methods  (4.1)  - 
(4.3)  when  applied  to  the  elliptic  and  parabolic  model  problems  (2.9) 
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and  (2.10)  with  the  k-line  and  k*k  block  structures  described  in 

i 

section  4.  We  limit  our  discussion  to  the  case  where  k  £  2;  although 
a  similar  argument  applies  when  k  =  1  (and  formulas  (6.4)  -  (6.9)  are 
valid  for  k  =  1),  we  have  not  in  that  instance  described  N.  We  begin 
by  showing  that  the  Jacobi  method  for  these  block  structures  satisfies 
the  assumptions  of  section  3.  After  p  is  determined  for  the  Jacobi 
method  it  is  easy  to  find  the  convergence  rates  of  the  Gauss-Seidel  and 
SOR  methods. 

Lemma  6.1.  Assumption  A1  holds  for  both  block  structures  and  both 
problems  if  h  is  sufficiently  small. 

Proof.  In  all  cases,  inspection  of  the  submatrices  (4.7)  of  A,  as 
given  by  (4.8),  shows  that  the  diagonal  elements  of  A  are  positive  and 
the  other  elements  are,  for  small  h,  nonpositive.  Therefore  N  is 
noonegative  and  A  and  M  are  M-matrices:  that  is, 

(6.1)  N  i  0,  M"1  Z  0,  and  A-1  2  0. 

Moreover,  A  is  irreducible  for  small  h.  A1  follows  from  Theorem  3.13 
of  [33] .  □ 


Lemma  6.2.  Assumption  A2  holds  for  both  block  structures  and  both 
problems  if  h  is  sufficiently  small. 

Proof.  This  follows  from  (6.1)  and  the  Perron-Frobenius  theory, 
see  Theorem  2.1  in  {33].  □ 
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We  remark  that  when  nonnegativity  of  M  *  or  A  *  fails,  A1  and  A2 
often  can  be  established  by  other  means,  for  example,  A1  holds  when  A 
is  positive  definite  and  N  is  nonnegative.  A2  follows  from  supposing 
that  M  is  positive  definite,  N  is  symmetric,  and  the  splitting 
satisfies  block,  property  A  (see  [1],  ( 33 1 ,  [35],  [25]),  for  then 
eigenvalues  of  M  are  real  and  occur  in  signed  pairs. 

Lemma  6.3.  Assumption  A3  holds  for  both  block  structures  and  both 
problems  if  h  is  sufficiently  small. 
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1.  1  T(i,s)  =  (l/k)(U,V)  +  E/2  +  hR/2, 

1  s 

IV 

with  E  satisfying  (5.7)  and 

IRI  ^  |e|J(l/k)(|U|,|V|)  +  IEI/2]. 

The  second  term  in  the  expansion  of  (NU,V)  is  appraised  in  the  same 
way,  to  yield 

(NU,V)  =  (2/k) (U,V)  +  E, 

1 E |  S  h(2+h|el  )[(|V U|,|V|)+(|U|,|V V|)J 
y  y 

+  h(|e|oo/k)(|U|,|V|)  +  h23k(l+h|e|oo/2)  [  (-AU,U)+(-AV,V)  ] . 
But  this  implies  the  inequality  of  A4,  with  D  given  by  (6.2),  □ 

Our  next  theorems  follow  immediately  from  these  lemmas  and 
Theorems  3.1  -  3.3. 

Theorem  6.5.  Let  p(kL)  and  p(kB)  denote  the  spectral  radii  for 
the  k-line  and  k*k  block  structures,  respectively,  of  the  block  Jacobi 
scheme  applied  to  the  elliptic  problem  (2.9).  Let  Aq  denote  the 
smallest  eigenvalue  of  the  problem 

t 

Lv  =  \v  in  Q,  v  =  0  on  3fl. 

Then 

p(kL)  =  1  -  (k/2)AQh2  +  o(h2), 

(6.4) 

p(kB)  =  1  -  (k/4)AQh2  +  o (h2 ) .  O 
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Theorem  6.6.  Let  p(SkL)  and  p(SkB)  denote  the  spectral  radii  for 
the  k-line  and  kxk  block  structures,  respectively,  of  the  block  Jacobi 
scheme  applied  to  the  parabolic  problem  (2.10),  and  suppose  that  PI 
holds.  Let 


A1  :=  min  (c(x,y)  :  (x,y)  €  Q} . 

/ 

Then 

p(SkL)  <  1  -  (k/2)A1ha  +  o(h“), 

(6.5) 

p(SkB)  S  1  -  (k/4)A1h“  +  o(ha), 
and  equality  holds  if  d  =  e  =  0,  so  that  SI  is  satisfied.  □ 


We  remark  that  the  character  "S"  is  to  remind  us  of  the  singular 
perturbation  nature  of  the  parabolic  equation. 

When  a  matrix  A  under  a  block  partition  satisfies  block  property 
A,  then  the  spectral  radii  p„c  of  the  Gauss-Seidel  method  (4.2)  and  p 

bo  U) 

of  the  SOR  method  (4.3)  are  determined  by  the  spectral  radius  p  of  the 
Jacobi  method  ([1],  (33,  chapter  4],  [35]): 

PGS  =  P2’  (PU)  +  W  -  1)2  =  w2P2Pw- 

Moreover,  p^  is  minimized  for  a  specific  u>: 

=  2/(1  +  (1  -  p2)1/2),  pb  =  -  1. 

With  the  block  structure  imposed  by  (4.5)  or  (4.6),  A  has  block 
property  A.  This  observation  proves  our  next  result. 
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2 

Corollary  6.7.  Let  A  represent  h  L^.  Then 
PGS(kL)  =  1  -  kAQh2  +  o(h2), 

(6.6) 

Pb(kL)  =  1  -  2(kAQ)1/2h  +  o(h), 
and 

PGS(kB)  =  1  -  (k/2)AQh2  +  o(h2), 

(6.7) 

Pb(kB)  =  1  -  (2kA0)1/2h  +  o(h). 

2 

Let  A  represent  h  1L.  Then 

PGS (SkL)  S  1  -  kAjh®  +  o(hff), 

(6.8) 

Pb(SkL)  <  1  -  2(kA1)1/2hCf/2  +  o(ha/2) , 

and 

PGS(SkB)  <  1  -  (k/2)A1ha  +  o(h“), 

(6.9) 

Pb(SkB)  %  1  -  (2kA1)1/2h0f/2  +  o(h“/2).  D 


7.  Other  operators.  In  this  section  we  extend  our  theory  to  cover  the 
more  general  operators  L  and  Z  defined  by 

Lu  :=  -  (au  )  -  (bu  )  +  du  +  eu  +  fu, 

xx  y  y  x  y 

(7.1) 

f 

£u  :=  cut  +  Lu. 

For  simplicity  we  have  excluded  terms  in  the  cross-derivative  u^y. 
Self-adjoint  operators  L  with  this  term  have  been  discussed  in  [23]. 

We  assume  for  convenience  that  a,  b,  c,  d,  e,  and  f  are  smooth 
functions  on  Q,  that  c  satisfies  (2.5),  and  that 


(7.2) 


a(x,y)  £  aQ  >  0,  b(x,y)  i  bQ  >  0,  f(x,y)  i  0  on  Q. 


L  is  uniformly  elliptic  by  the  strict  positivity  of  a  and  b,  and 
satisfies  a  maximum  principle  by  virtue  of  the  nonnegativity  of  f. 

As  in  section  2,  we  let  U  be  a  mesh  vector  on  the  mesh  points 
defined  by  (2.3).  At  points  (x^y^)  of  we  define 


a . 


i+^,j 


■•=  (a. 


i.J 


+  a 


i+l.j 


J/2,  b. 


:=  (b 


+  b. 


i.j+1 


)/2. 


The  discrete  approximations  to  (7.1)  are  then 


(7.3) 


and 


‘V’i.j  -  '  '  ui,j>  '  '  ui-U)1/h 

*  -  “i.j’  -  bi.j->,(ui,j  '  "i,j-.»/b2 

+  ‘ij'Vl.J  -  «i-idW<W  * 


(7.4)  IVli.j  -  *  'ViJ' 


It  is  not  difficult  to  see  that  the  machinery  of  section  3  still 
works.  The  main  theorem  of  124],  which  relates  the  minimal  eigenvalues 
of  (3.1)  and  (3.5),  is  easy  to  establish  with  L  and  Lfa  given  by  (7.1) 
and  (7.3),  respectively.  Consequently  Theorems  3.1,  3.2,  and  3.3 
apply,  mutatis  mutandis,  to  splittings  of  the  matrix  A  arising  from 
(7.3)  or  (7.4). 

Now  we  must  determine  q  for  the  Jacobi  scheme,  using  either  of  the 
block  structures  (4.5)  or  (4.6), 

For  the  k-line  structure,  a  direct  computation  yields 


(KU,V)  =  I  Ib±>ks+!i 


♦  I  [b 


i  ,ks+*j 


+  he 

-  he 


i,ks+l^^i,ksVi,ks+l 
i,ks^^i,ks+l^i,ks  ’ 
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where  the  sum  is  over  lSiSP,  0£s£  Q-l.  Consider  a  term 


T(i,s)  (bi>ks+^  +  hei>ks+1/2lui>ksVi>ks+l 

~  ^i,ks+j^i,ks^i,ks+l  +  ^^ei ,ks+l^^i ,ks^i ,ks+l 
+  ^i,ks+Jj  **i,ks+j  ^i,ks^i,ks+l* 

The  factor  in  square  brackets  in  the  last  term  above  is  bounded  by 
hklVybl^,  because  b  is  smooth.  Proceeding  as  in  the  proof  of  Lemma 
6.4,  we  establish  the  validity  of  A4  with 


(7.5)  q  =  2b/k,  D  =  max  Uklbl.+OCh) ,  2|V vblw+|e|0#/k}. 


Observe  that  the  variable  coefficient  b  has  led  to  a  variable  q. 
In  the  same  way,  for  the  kxk  block  scheme  we  obtain 


q  *  (2a  +  2b)/k, 

(7.6) 

D  =  max  {3k|a |  +3k|b|  +0(h) ,  2|V  a|„+2|V b|a>+(|d|0  +  |e|a>)/k} . 


We  collect  our  results  in  the  following  two  theorems. 


Theorem  7.1.  Let  p(kL)  and  p(kB)  denote  the  spectral  radii  for 
the  horizontal  k-line  and  k*k  block  structures,  respectively,  of  the 
block  Jacobi  scheme  applied  to  the  elliptic  problem  (2.9)  with  L^  given 
by  (7.3).  Let  denote  the  smallest  eigenvalue  of  the  problem 

f 

Lv  =  ybv  in  0,  v  =  0  on  30, 

and  let  Tq  denote  the  smallest  eigenvalue  of  the  problem 
Lv  =  y(a+b)v  in  0,  v  =  0  on  30. 
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Then 


p(kL)  =  1  -  (k/2)fQh2  +  o(h2), 

(7.7) 

p(kB)  =  1  -  (k/2)rQh2  +  o(h2).  0 

Theorem  7.2.  Let  p(SkL)  and  p(SkB)  denote  the  spectral  radii  for 
the  horizontal  k-line  and  k*k  block  structures,  respectively,  of  the 
block  Jacobi  scheme  applied  to  the  parabolic  problem  (2.10)  with 
given  by  (7.4),  and  suppose  that  PI  holds.  Let 

fj  :=  min  {c(x,y)/b(x,y)  :  (x,y)  e  fi}, 

Tj  :=  min  {c(x,y)/(a(x,y)  +  b(x,y))  :  (x,y)  €  fl} . 

Then 

p(SkL)  S  1  -  (k/2)f1h°  +  o(h°) , 

(7.8) 

p(SkB)  5  1  -  (k/2)rih“  +  o(ha), 
and  equality  holds  if  d  =  e  =  0,  so  that  SI  is  satisfied.  □ 

The  nonzero  pattern  of  A  is  the  same  whether  A  arises  from  (2.6) 
or  (7.3).  In  the  more  general  case,  then,  A  retains  block  property  A 
for  both  the  k-line  and  kxk  block  partitions.  Consequently  the 
analogue  of  Corollary  6.7  is  valid.  We  leave  a  statement  of  this 
Corollary  7.3  to  the  reader. 
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8.  Other  domains.  Extension  of  our  results  from  two  to  m  space 

* 

dimensions  is  straightforward.  We  sketch  this  for  the  model  problems 
set  in  the  m-dimensional  unit  cube,  and  then  return  to  the 
two-dimensional  setting  to  discuss  domains  other  than  the  unit  square. 
Treatment  of  general  domains  in  higher  dimensions  is  similar,  while 
more  general  operators  on  these  domains  can  be  handled  as  outlined  in 
section  7. 

Impose  a  uniform  mesh  of  si2e  h  =  1/(P+1)  on  the  unit  cube 

Q(m)  :=  {x  =  (x. ,  . . ,x  )  €  Rm  :  0  <  x.  <  1}. 

i  m  ~  x 

Let  p  be  a  multi-index  p  =  (pj,  . .,Pm)  €  Z®,  and  let  y(i)  be  the 
multi-index  whose  ith  component  is  1  and  whose  other  components  are  0. 
By  Xp  we  mean  the  point  x  =  (p^h,  ..jp^h)  in  Rm.  Hence  putting 

B  :=  {p  €  Z®  :  1  £  P^  £  P},  B  :=  {p  e  Z®  :  0  S  p.  <  P+l} 
allows  us  to  write 

OjjOn)  =  {Xp  :  p  €  B}, 

3fl^(m)  =  {xp  :  p  €  B  and  at  least  one  p.  =  0  or  =  P+lj. 

2  2 

With  the  m-dimensional  Laplacian  A(m)  :  =  1^  9  /9x^  ,  our  model 
operators  are 

4 

Lu  :=  -A(m)u  +  lT_j  d^9u/9x^, 

(8.1) 

£u  :=  c9u/9t  +  Lu; 

we  suppose  that  c  and  all  d^  are  smooth  and  that  c  satisfies  (2.5). 
Discretization  of  these  operators  is  done  as  in  (2.6)  and  (2.7).  Let 
U  =  (Up)  be  a  mesh  vector.  The  approximation  of  A(m)  is  given  by 
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(8.2) 


and  the 

(8.3) 

(8.4) 


iV-Mp :=  zi  (VyU)  ■  2Ue  *  Vv(i))/h2' 

* 

discrete  operators  corresponding  to  (8.1)  are 

iV'p  *  zi  "i,p(VY(i)  •  Vv(i))/(2h)' 

IV]P ;=  (cp/T)ue  *  <W 


In  m  dimensions  there  are  m  obvious  block  partitions.  Suppose  for 
example  that  m  =  3.  In  the  k-plane  block  structure,  each  block  of  U 
comprises  the  unknowns  Up  associated  with  the  points  Xp  on  k 
consecutive  planes.  Indexing  the  blocks  by  s,  we  have 


Us  :=  :  p  €  B,  k(s-l)  <  P3  S  ks}. 

In  like  fashion,  for  blocks  of  k*k  lines  the  basic  subblock  is 

Ur  s  :=  {Up  :  p  €  B,  k(r-l)  <  P2  g  kr,  k(s-l)  <  P3  £  ks}, 

and  kxk*k  blocks  are  given  by 

0  .  ,  :=  {UR  :  P  €  B,  k(r-l)  <  p.  £  kr,  k(s-l)  <  P,  £  ks, 
r,a,t  p  1  * 

k(t-l)  <  p3  <  kt}. 

Let  us  agree  to  call  a  basic  block  an  s-slice  of  U  if  we  partition  U  by 
imposing  restrictions  of  the  form  k(s^-l)  <  P^  £  ks^  on  the  indices 
^■-s+I’  ^m*  b^is  notation,  it  is  easy  to  state  and  prove  the 

■-dimensional  version  of  Theorem  4.1. 


Theorem  8.1.  Let  k  £  2.  Denote  by  the  operator  acting  as  N  in 

2  2 

the  Xj-direction.  Let  A  represent  h  L^  or  h  £^.  For  the  block  Jacobi 
scheme  (4.1)  based  on  s-slice  blocks, 

» -  z;™—i  »j  *  °<h>-  ° 
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Observe  that  s-slice  decompositions  preserve  block  property  A. 
Assertions  similar  to  Lemmas  6.1  -  6.4  are  readily  demonstrated;  the 
following  lemma  collects  the  results. 

2  2 

Lemma  8.2.  Let  k  i  2  and  let  A  correspond  to  h  L^  or  h 
Assumptions  A1  -  A4  are  satisfied  by  the  block  Jacobi  method  (4.1) 
based  on  s-slices.  Specifically,  for  1  £  s  S  m, 

INI k  3  NQ  =  2s/k  +  0(h)  S  s  +  1,  q  =  2s/k, 

D  =  max  l3sk+0(h),  ^_m_s+1  Id^l^k}, 

(8.5) 

|EI  S  hD[L”=1  (IV  U|,|V|)  +  l”=1  (lUlJVjVl)  +  (|U|,|V|)] 

+  h^K-^WU.D)  +  (-^WV.V)].  O 

Now  the  machinery  of  section  3  grinds  out  theorems  like  those  of 

section  6.  Rather  than  turn  the  crank,  we  choose  to  consider  problems 

2 

with  the  model  operators  (2.1)  set  in  more  general  domains  0  of  R  . 

We  begin  by  describing  ft  and  Assume  that  Q  is  a  bounded 

domain  in  R  with  Lipschitz  boundary  and  that  locally  0  lies  always 
on  one  side  of  8(1.  This  last  condition  ensures  that  Q  has  no  internal 
cusps.  Boundedness  implies  that  (1  has  "leftmost"  and  "bottommost" 

—  4  m 

tangents  x  =  xQ  :=  min  {x  :  (x,y)  €  fl},  y  =  yQ  :=  min  {y  :  (x,y)  e  0}. 

2 

Choose  h  >  0  and  impose  on  R  a  grid  of  lines 

(8.6)  x  =  x .  :=  xQ  +  ih,  y  =  :  =  yQ  +  jh; 

the  intersections  (x^.y^)  are  called  grid  points ■  Define  to  be  the 
collection  of  all  the  grid  points  (x^y^)  6  Q  together  with  all  the 


33 


points  of  intersection  of  the  grid  lines  (8.6)  with  3Q.  Then  80^  := 
n  80,  and  0^  consists  of  those  points  of  0^  that  lie  in  0.  The 
points  of  0^  are  called  mesh  points. 

The  points  of  0^  are  conveniently  viewed  as  being  of  two  types: 
the  four  nearest  neighbors  of  a  regular  point  are  themselves  grid 
points  (x^.y^),  and  the  other  points  of  0^  are  called  irregular.  We 
consider  finite  difference  approximations  of  the  operators  (2.1)  that 
differ  from  the  earlier  constructions  (2.6)  and  (2.7)  only  at  irregular 
points.  Any  of  a  number  of  approximations  at  an  irregular  point  will 
suffice  for  our  purposes;  it  is  only  necessary  that  the  approximation 
be  at  least  an  interpolation  of  degree  0  (see  [13,  pp.  199  ff.],  [6]). 
Proceeding  as  before  brings  us  to  a  linear  system  (2.11),  to  which  the 
block  iterative  methods  of  section  4  can  be  applied. 

To  use  the  theory  of  section  3,  we  need  to  establish  the  relations 

(NU,V)  =  (qU.V)  +  E, 

(8.7)  |E|  S  hD[(|V  U|+|V  U|,|V|)  +  (|U| , |V  V|+|V  V|)  +  (|U|,|V|)] 

*  y  a  y 

+  h^K-AjU.U)  +  (-AjV.V)] 

of  A4  for  mesh  vectors  U,  V  that  vanish  on  80^.  It  is  clear  from  the 
argument  of  Lemma  6.4  that  (8.7)  remains  true  if  Q  is  composed  of  "grid 

4 

rectangles,"  whose  sides  are  grid  lines  (8.6).  All  that  remains  is  to 
prove  (8.7)  for  more  general  domains.  But  this  will  follow  if  we  can 
show  that  Q  is  almost  a  union  of  grid  rectangles.  To  this  end,  note 
that  Q,  being  bounded,  is  contained  in  some  rectangle 

R  :=  [(x,y)  e  R2  :  xQ  S  x  S  xp+1,  yQ  S  V  *  yp+1). 
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) 


Pick  an  integer  k  £  2  and  subdivide  R  into  closed  subrectangles 


fl  is  the  union  of  rectangles  R  that  lie  entirely  within  Q,  and 

r ,  s 

fragments  of  such  rectangles.  Let  R  be  the  interior  of  the  set 


and  for  any  subset  G  of  R  let  (U,V)g  denote  the  inner  product  of  mesh 
vectors  U  and  V  over  the  mesh  points  in  G.  Because  N  is  bounded,  (8.7) 
is  an  easy  consequence  of  the  following  lemma,  whose  proof  is  similar 
to  the  proof  of  inequality  (13)  in  [11]. 


Lemma  8.3.  Let  U  and  V  be  mesh  vectors  that  vanish  on  3(1^.  Then 
(8.8)  !(U,V)fl  -  (U,V)RJ  £  2h2w2 | (-A^U ,U)  +  (-A,V,V)] 
whenever  h  %  hQ,  for  some  constants  w  and  that  depend  only  on  ft. 


Proof.  Divide  8fl  into  a  finite  number  of  pieces  for  which  the 
angle  of  the  tangent  with  either  the  x-  or  y-axis  exceeds  some  positive 
value  (say  30°).  For  instance,  let  B  be  a  piece  of  the  boundary  that 
is  this  steep  with  respect  to  the  x-axis.  Let  S  be  the  horizontal 

0  4 

strip  of  Q  that  abuts  B  and  is  k  grid  lines  high,  so  that  S  is  for  some 

fixed  s  the  union  of  rectangles  R  and  at  most  two  pieces  composed  of 

r ,  s 

fragments  of  such  rectangles.  Denote  by  F  the  piece  touching  B  —  say 
the  leftmost  piece.  The  smoothness  of  9ft  ensures  that  there  are 
positive  constants  h^  and  w,  depending  only  on  ft,  so  that  the  x-width 
of  F  is  at  most  wh  when  h  $  h^.  See  Figure  l. 
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Now  consider  a  horizontal  mesh  line  y  =  y^  in  this  strip  S;  let 

be  the  leftmost  and  the  rightmost  mesh  points  on  the  line,  and 

let  x^  be  the  first  point  in  the  leftmost  subrectangle  in  S.  Observe 

that  |xj  -  xfl|  £  wh.  Because  U  vanishes  on  B,  at  each  point  xi  between 

x  and  x,  we  have 
a  f 

Ui,j  "  *0=3  ^Uo+l,j  "  Uo,j>> 

and  similarly  for  V.  Hence  |U.  .|  i  hi*  *  |V  U  . |,  so 

i,j  o=a  x  o,j  ’ 

S"2<£1  'Vo,jU 

$  h2w(l*~}  |V  U_  ■|2)1/2(I*I1  |V  V  .|2)1/2 
o=a  x  a,j  *  v  o=a  x  ct,j 

5  t^wl*!1  [  (V  U  .|2  +  |V  V  .  | 2 ] /2 

0=a  x  o,j 1  x  o,j 

6  h2wl^1  [|V  U  -|2  +  |V  V  .  |2]/2. 

o= a  x  <j,j 1  x  cr, j 

Now  sum  over  x.  between  x  and  xr  to  get 
i  a  f 


|Iiui.jvijl  SIi 

s  h2“2^I  I'Voj'2  +  'Vo,/"2- 

and  sum  this  last  inequality  over  j  to  see  that 

IOJ,v)Ft  s  tV[ivxui^s  ♦  iV'h,s1/2- 

Repeat  this  over  every  fragment  F  and  boundary  piece  B.  Then  each 

subrectangle  R  within  f)  is  covered  at  most  four  times,  so 
r  ,s 

|(U,V)fl  -  (U,V)R|  S  2h2w2{lvu#2  +  ivulj  +  IVXVI2  +  lv  vl2]; 


(8.8)  follows  from  Lemma  2.1.  □ 
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9.  How  general  is  the  method?  We  expect  that  oar  techniques  will 
permit  us  to  analyze  any  block  iterative  scheme  in  which  the  blocks 
have  a  regular  pattern,  so  that  NU  constitutes  an  orderly,  weighted 
sampling  of  any  mesh  vector  U.  Evidently  most  finite  difference 
approximations  on  nice  meshes  give  rise  to  such  operators  N. 

For  instance,  on  uniform  meshes  in  Rm  it  will  be  true  that 

h2[Lbu|p  *  zj  VVV 

where  £  runs  over  the  set  of  indices  E((3)  of  "neighboring"  mesh  points 
of  the  point  x^.  Each  smooth  enough  coefficient  Ag(x^)  will  satisfy 

A 

Ap(x^)  =  A^(Xp)  +  0(h) 

for  some  smooth  function  A^,  and  will  be  a  linear  combination  of  the 

coefficients  of  the  differential  operator  L  to  which  is  an 

approximation.  Because  N  is  derived  from  the  matrix  A  representing 
2 

h  L^,  we  will  have 

(9.1)  fNU]p  =  2^  Bp(x^)U^. 

If  N  is  sufficiently  regular,  then  there  will  be  some  regularly 
distributed  subset  S  of  points  of  0^  so  that 

0  xp*S 

(9.2)  Bp(x^)  =  ^ 

A 

B^(xp)  +  0(h)  |  e  Z(0),  Xp  e  S. 

Consider  now  a  typical  term  in  (NU,V).  At  a  point  Xp  of  S, 

irolpVp  =  it  BpUjWjVp. 
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It  follows  from  (9.1)  and  (9.2)  that  if  U  and  V  are  smooth  then 


|mVe s  h  VVW 

hence 

(9.3)  (NU,V)  S  (qU,V), 


where  q  accounts  for  the  terms  and  the  relative  cardinalities  of 
Z(0),  S,  and  The  approximate  inequality  of  (9.3)  indicates  that  N 

is  a  weak  multiplication  operator.  Development  of  an  estimate  of  the 
error  in  (9.3)  proceeds  as  in  section  5  and  in  the  proof  of  Lemma  6.4. 

In  effect,  our  view  of  N  as  a  weak  multiplication  operator  has 
already  been  used  in  [23,  section  7],  where  a  splitting  somewhat 
different  from  the  usual  k-line  block  Jacobi  scheme  is  treated. 


Rectangular  meshes,  which  have  uniform  spacing  h^  in  the  x^direction, 
can  also  be  handled,  as  can  other  boundary  conditions.  We  see  then 
that  this  viewpoint  unifies  the  derivation  of  the  convergence  rates  of 
block  iterative  methods  for  elliptic  and  parabolic  finite  difference 
equations.  In  fact,  the  technique  will  also  yield  estimates  of  the 
rates  of  convergence  of  block  relaxation  methods  applied  to  the 
matrices  arising  from  finite  element  approximations.  But  finite 
elements  are  powerful  in  part  because  they  admit  irregular  partitions 
of  Q.  For  such  partitions  it  is  not  apparent  how  to  group  the,  unknowns 
so  that  a  systematic  block  iterative  scheme  is  easy  to  implement.  We 
direct  the  reader  to  [4]  for  an  example  of  a  successful  application  of 
these  ideas  in  a  simple  case,  where  the  iterative  method  was  easy  to 


program. 
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Figure  1.  The  strip  S  of  Q. 

F  is  the  union  of  the  two  leftmost  fragments. 


43 


